"""
Phase 10: Monetary Union Comparison — Beyond the Eurozone
==========================================================
Tests whether the eurozone amplification result generalizes to other
monetary unions: CFA Franc Zone (WAEMU + CEMAC), Eastern Caribbean
Currency Union (ECCU), and Common Monetary Area (CMA).

Key hypothesis: if EMU amplifies aging-driven imbalances and CFA amplifies
youth-driven imbalances, the mechanism is about removing the exchange rate
absorber generically — a much stronger policy claim.

Sections:
1. Within-union CA regressions (each union separately)
2. Pooled monetary union analysis (union dummies + Z×union interactions)
3. Within-union Z deviations → CA deviations
4. CFA zone deep dive (anchor country analysis, projections)
5. Cross-union comparison summary
6. Regime strain projections for non-EMU unions
"""

import pandas as pd
import numpy as np
from pathlib import Path
import sys
import warnings
warnings.filterwarnings('ignore')

PROJECT_DIR = Path(__file__).resolve().parent.parent
ROOT_DIR = PROJECT_DIR.parent
sys.path.insert(0, str(ROOT_DIR / "multilateral" / "src"))
from model import PanelGLS

DATA = PROJECT_DIR / "data" / "processed"
OUT_TABLES = PROJECT_DIR / "output" / "tables"
OUT_TABLES.mkdir(parents=True, exist_ok=True)

MULTILATERAL_DATA = ROOT_DIR / "multilateral" / "followup" / "data" / "processed"

# ── Monetary Union Definitions ────────────────────────────────────────

# Eurozone (for comparison baseline)
EUROZONE_JOIN = {
    'AUT': 1999, 'BEL': 1999, 'FIN': 1999, 'FRA': 1999, 'DEU': 1999,
    'IRL': 1999, 'ITA': 1999, 'LUX': 1999, 'NLD': 1999, 'PRT': 1999,
    'ESP': 1999, 'GRC': 2001, 'SVN': 2007, 'CYP': 2008, 'MLT': 2008,
    'SVK': 2009, 'EST': 2011, 'LVA': 2014, 'LTU': 2015,
}

# CFA Franc Zone — WAEMU (West African Economic and Monetary Union)
# Pegged to French franc (1945), then euro (1999). Anchor: FRA/ECB
# Note: GNB joined in 1997
WAEMU_JOIN = {
    'BEN': 1960, 'BFA': 1960, 'CIV': 1960, 'MLI': 1984,  # Mali rejoined 1984
    'NER': 1960, 'SEN': 1960, 'TGO': 1960, 'GNB': 1997,
}

# CFA Franc Zone — CEMAC (Central African Economic and Monetary Community)
CEMAC_JOIN = {
    'CMR': 1960, 'CAF': 1960, 'TCD': 1960, 'COG': 1960,
    'GNQ': 1985, 'GAB': 1960,
}

# Combined CFA
CFA_JOIN = {**WAEMU_JOIN, **CEMAC_JOIN}
CFA_ISO3 = set(CFA_JOIN.keys())
WAEMU_ISO3 = set(WAEMU_JOIN.keys())
CEMAC_ISO3 = set(CEMAC_JOIN.keys())

# Eastern Caribbean Currency Union (ECCU)
# EC dollar pegged to USD since 1976. Anchor: USA
ECCU_JOIN = {
    'ATG': 1965, 'DMA': 1965, 'GRD': 1965, 'KNA': 1965,
    'LCA': 1965, 'VCT': 1965,
    # AIA and MSR have limited data; include if available
    'AIA': 1965, 'MSR': 1965,
}
ECCU_ISO3 = set(ECCU_JOIN.keys())

# Common Monetary Area (CMA)
# Southern Africa: rand-linked currencies. Anchor: ZAF
CMA_JOIN = {
    'ZAF': 1974, 'LSO': 1974, 'SWZ': 1974, 'NAM': 1990,
}
CMA_ISO3 = set(CMA_JOIN.keys())

EUROZONE_ISO3 = set(EUROZONE_JOIN.keys())

# All unions with metadata
UNIONS = {
    'EMU': {
        'members': EUROZONE_JOIN, 'iso3': EUROZONE_ISO3,
        'anchor': 'DEU', 'anchor_currency': 'EUR',
        'label': 'Eurozone', 'peg_target': 'euro',
    },
    'CFA': {
        'members': CFA_JOIN, 'iso3': CFA_ISO3,
        'anchor': 'FRA', 'anchor_currency': 'EUR (via CFA franc)',
        'label': 'CFA Franc Zone', 'peg_target': 'euro',
        'subgroups': {'WAEMU': WAEMU_ISO3, 'CEMAC': CEMAC_ISO3},
    },
    'ECCU': {
        'members': ECCU_JOIN, 'iso3': ECCU_ISO3,
        'anchor': 'USA', 'anchor_currency': 'USD (via EC dollar)',
        'label': 'Eastern Caribbean (ECCU)', 'peg_target': 'USD',
    },
    'CMA': {
        'members': CMA_JOIN, 'iso3': CMA_ISO3,
        'anchor': 'ZAF', 'anchor_currency': 'ZAR',
        'label': 'Common Monetary Area', 'peg_target': 'ZAR',
    },
}

OECD = {
    'AUS', 'AUT', 'BEL', 'CAN', 'CHL', 'COL', 'CRI', 'CZE', 'DNK', 'EST',
    'FIN', 'FRA', 'DEU', 'GRC', 'HUN', 'ISL', 'IRL', 'ISR', 'ITA', 'JPN',
    'KOR', 'LVA', 'LTU', 'LUX', 'MEX', 'NLD', 'NZL', 'NOR', 'POL', 'PRT',
    'SVK', 'SVN', 'ESP', 'SWE', 'CHE', 'TUR', 'GBR', 'USA',
}


def stars(p):
    if p < 0.01: return '***'
    if p < 0.05: return '**'
    if p < 0.1: return '*'
    return ''


def fmt(val, se, p):
    s = stars(p)
    return f"{val:.4f}{s}", f"({se:.4f})"


def run_panel_gls(df, y_var, x_vars, label):
    """Run PanelGLS and return results dict."""
    cols = [y_var] + x_vars + ['iso3', 'year']
    sub = df[cols].dropna()
    if len(sub) < 30:
        print(f"  {label}: insufficient obs ({len(sub)}), skipping")
        return None

    n_countries = sub['iso3'].nunique()
    if n_countries < 3:
        print(f"  {label}: insufficient countries ({n_countries}), skipping")
        return None

    gls = PanelGLS()
    y = sub[y_var].values
    X = sub[x_vars].values
    try:
        gls.fit(y, X, sub['iso3'].values, sub['year'].values)
    except Exception as e:
        print(f"  {label}: GLS failed ({e}), skipping")
        return None

    result = {
        'model': label,
        'n_obs': gls.n_obs,
        'n_countries': gls.n_countries,
        'r_squared': gls.r_squared,
        'rho': gls.rho,
    }
    print(f"\n  {label} (N={gls.n_obs}, {gls.n_countries} countries, "
          f"R²={gls.r_squared:.4f})")
    for i, name in enumerate(x_vars):
        sig = stars(gls.pvalues[i])
        print(f"    {name:30s} {gls.beta[i]:8.4f} ({gls.se[i]:.4f}) {sig}")
        result[f'{name}_coef'] = gls.beta[i]
        result[f'{name}_se'] = gls.se[i]
        result[f'{name}_p'] = gls.pvalues[i]

    result['_gls'] = gls
    result['_x_vars'] = x_vars
    return result


def write_table(results, filename, title, note=None):
    """Write regression results as markdown table."""
    if not results:
        return

    lines = [f"# {title}\n"]

    all_vars = []
    for r in results:
        for k in r:
            if k.endswith('_coef'):
                vname = k.replace('_coef', '')
                if vname not in all_vars:
                    all_vars.append(vname)

    model_labels = [r['model'] for r in results]
    header = "| Variable | " + " | ".join(model_labels) + " |"
    sep = "|:---|" + "|".join(["---:" for _ in results]) + "|"
    lines.append(header)
    lines.append(sep)

    for var in all_vars:
        coef_row = f"| {var} |"
        se_row = "| |"
        for r in results:
            if f'{var}_coef' in r:
                c, s = fmt(r[f'{var}_coef'], r[f'{var}_se'], r[f'{var}_p'])
                coef_row += f" {c} |"
                se_row += f" {s} |"
            else:
                coef_row += " |"
                se_row += " |"
        lines.append(coef_row)
        lines.append(se_row)

    lines.append("|:---|" + "|".join(["---:" for _ in results]) + "|")
    n_row = "| N |"
    r2_row = "| R² |"
    nc_row = "| Countries |"
    for r in results:
        n_row += f" {r['n_obs']} |"
        r2_row += f" {r['r_squared']:.4f} |"
        nc_row += f" {r['n_countries']} |"
    lines.append(n_row)
    lines.append(r2_row)
    lines.append(nc_row)

    if note:
        lines.append(f"\n{note}")
    else:
        lines.append("\n*Panel GLS with country and year fixed effects. "
                     "Standard errors in parentheses.*")
        lines.append("*\\*p<0.1, \\*\\*p<0.05, \\*\\*\\*p<0.01*")

    path = OUT_TABLES / filename
    path.write_text('\n'.join(lines))
    print(f"\n  Saved: {path}")


def get_union_df(df, union_members):
    """Filter to union members post-join-year."""
    rows = []
    for iso3, join_yr in union_members.items():
        mask = (df['iso3'] == iso3) & (df['year'] >= join_yr)
        rows.append(df[mask])
    if not rows:
        return pd.DataFrame()
    return pd.concat(rows, ignore_index=True)


# ── 1. Within-Union CA Regressions ───────────────────────────────────

def within_union_ca(df):
    """Run Z → CA regressions within each monetary union."""
    print("\n" + "=" * 60)
    print("1. WITHIN-UNION CA REGRESSIONS")
    print("=" * 60)

    z_vars = ['Z_1', 'Z_2', 'Z_3']
    controls = ['fiscal_bal_gdp', 'nfa_gdp_lag', 'rgdp_growth']
    all_results = {}

    for union_key, union_info in UNIONS.items():
        print(f"\n  --- {union_info['label']} ---")
        udf = get_union_df(df, union_info['members'])
        if len(udf) == 0:
            print(f"  No data for {union_key}")
            continue

        n_ca = udf['ca_gdp'].notna().sum()
        print(f"  Sample: {len(udf)} obs, {udf['iso3'].nunique()} countries, "
              f"{n_ca} with CA data")

        results = []

        # M1: Z levels → CA (full controls)
        r = run_panel_gls(udf, 'ca_gdp', z_vars + controls,
                          f'{union_key}: Z+ctrl')
        if r:
            results.append(r)

        # M2: Z only (no controls) — useful for small samples
        r = run_panel_gls(udf, 'ca_gdp', z_vars,
                          f'{union_key}: Z only')
        if r:
            results.append(r)

        # M3: Age decomposition
        age_vars = ['old_dep', 'youth_dep']
        avail_age = [v for v in age_vars if v in udf.columns]
        if avail_age:
            r = run_panel_gls(udf, 'ca_gdp', avail_age + controls,
                              f'{union_key}: Age')
            if r:
                results.append(r)

        all_results[union_key] = results

    # Write combined comparison table (Z+ctrl results across unions)
    comparison = []
    for union_key in ['EMU', 'CFA', 'ECCU', 'CMA']:
        if union_key in all_results:
            for r in all_results[union_key]:
                if r and 'Z+ctrl' in r['model']:
                    comparison.append(r)
                    break
            else:
                # Fall back to Z-only if Z+ctrl failed
                for r in all_results[union_key]:
                    if r and 'Z only' in r['model']:
                        comparison.append(r)
                        break

    if comparison:
        write_table(comparison, "phase10_within_union_ca.md",
                    "Within-Union CA Regressions: Cross-Union Comparison",
                    note=("*Panel GLS with country and year fixed effects. "
                          "Each union estimated separately, members post-accession only. "
                          "Controls: fiscal_bal_gdp, nfa_gdp_lag, rgdp_growth.*\n"
                          "*\\*p<0.1, \\*\\*p<0.05, \\*\\*\\*p<0.01*"))

    # Also write individual union tables with all specs
    for union_key, results in all_results.items():
        if results:
            write_table(results,
                        f"phase10_{union_key.lower()}_ca_detail.md",
                        f"{UNIONS[union_key]['label']}: CA Regression Detail")

    return all_results


# ── 2. Pooled Monetary Union Analysis ────────────────────────────────

def pooled_union_analysis(df):
    """Pooled regression with union dummies and Z×union interactions."""
    print("\n" + "=" * 60)
    print("2. POOLED MONETARY UNION ANALYSIS")
    print("=" * 60)

    z_vars = ['Z_1', 'Z_2', 'Z_3']
    controls = ['fiscal_bal_gdp', 'nfa_gdp_lag', 'rgdp_growth']

    # Create union membership indicators (post-join)
    pdf = df.copy()
    for union_key, union_info in UNIONS.items():
        col = f'in_{union_key.lower()}'
        pdf[col] = 0.0
        for iso3, join_yr in union_info['members'].items():
            mask = (pdf['iso3'] == iso3) & (pdf['year'] >= join_yr)
            pdf.loc[mask, col] = 1.0

    # Combined "any monetary union" dummy
    union_cols = [f'in_{k.lower()}' for k in UNIONS]
    pdf['in_any_union'] = (pdf[union_cols].sum(axis=1) > 0).astype(float)

    results = []

    # M1: Baseline — Z + controls (full sample)
    r = run_panel_gls(pdf, 'ca_gdp', z_vars + controls, 'Baseline')
    if r:
        results.append(r)

    # M2: Add "any union" dummy
    r = run_panel_gls(pdf, 'ca_gdp', z_vars + controls + ['in_any_union'],
                      '+ Union dummy')
    if r:
        results.append(r)

    # M3: Z×union interactions (any union)
    for z in z_vars:
        pdf[f'{z}_x_union'] = pdf[z] * pdf['in_any_union']
    z_x_union = [f'{z}_x_union' for z in z_vars]
    r = run_panel_gls(pdf, 'ca_gdp',
                      z_vars + controls + ['in_any_union'] + z_x_union,
                      'Z × Union')
    if r:
        results.append(r)

    # M4: Separate union dummies
    avail_union_cols = [c for c in union_cols if pdf[c].sum() > 0]
    r = run_panel_gls(pdf, 'ca_gdp',
                      z_vars + controls + avail_union_cols,
                      'Union dummies')
    if r:
        results.append(r)

    # M5: Z₁ × each union interaction
    for union_key in UNIONS:
        col = f'in_{union_key.lower()}'
        if pdf[col].sum() > 0:
            pdf[f'Z_1_x_{union_key.lower()}'] = pdf['Z_1'] * pdf[col]
    z1_x_each = [f'Z_1_x_{k.lower()}' for k in UNIONS
                 if f'Z_1_x_{k.lower()}' in pdf.columns]
    r = run_panel_gls(pdf, 'ca_gdp',
                      z_vars + controls + avail_union_cols + z1_x_each,
                      'Z₁ × each union')
    if r:
        results.append(r)

    write_table(results, "phase10_pooled_unions.md",
                "Pooled Monetary Union Analysis: Z → CA with Union Interactions",
                note=("*Panel GLS with country and year fixed effects. "
                      "Full sample. Union membership = 1 for member-years "
                      "post-accession. Z×Union tests whether monetary union "
                      "membership amplifies the demographic channel.*\n"
                      "*\\*p<0.1, \\*\\*p<0.05, \\*\\*\\*p<0.01*"))

    return pdf, results


# ── 3. Within-Union Z Deviations → CA Deviations ─────────────────────

def within_union_deviations(df):
    """Compute Z deviations from union mean and test deviation → CA deviation."""
    print("\n" + "=" * 60)
    print("3. WITHIN-UNION Z DEVIATIONS → CA DEVIATIONS")
    print("=" * 60)

    z_vars = ['Z_1', 'Z_2', 'Z_3']
    controls = ['fiscal_bal_gdp', 'nfa_gdp_lag', 'rgdp_growth']
    all_results = []

    for union_key, union_info in UNIONS.items():
        print(f"\n  --- {union_info['label']} ---")
        udf = get_union_df(df, union_info['members'])
        if len(udf) < 30:
            print(f"  Insufficient data ({len(udf)} obs)")
            continue

        # Compute union-year means
        u_means = udf.groupby('year')[z_vars + ['ca_gdp']].mean()
        u_means.columns = [f'{c}_union_mean' for c in u_means.columns]
        udf = udf.merge(u_means, on='year', how='left')

        for z in z_vars:
            udf[f'{z}_dev'] = udf[z] - udf[f'{z}_union_mean']
        udf['ca_dev'] = udf['ca_gdp'] - udf['ca_gdp_union_mean']

        # Distance from anchor
        anchor = union_info['anchor']
        if anchor in udf['iso3'].values:
            anchor_z = udf[udf['iso3'] == anchor][['year'] + z_vars].set_index('year')
            anchor_z.columns = [f'{c}_anchor' for c in anchor_z.columns]
            udf = udf.merge(anchor_z, on='year', how='left')
            for z in z_vars:
                udf[f'{z}_dist_anchor'] = udf[z] - udf[f'{z}_anchor']
        else:
            print(f"  Anchor {anchor} not in union data")

        # Z deviation → CA deviation
        z_dev_vars = [f'{z}_dev' for z in z_vars]
        r = run_panel_gls(udf, 'ca_dev', z_dev_vars,
                          f'{union_key}: Z dev')
        if r:
            all_results.append(r)

        # Z deviation → CA deviation + controls
        r = run_panel_gls(udf, 'ca_dev', z_dev_vars + controls,
                          f'{union_key}: Z dev+ctrl')
        if r:
            all_results.append(r)

        # Distance from anchor → CA (excl anchor)
        z_dist_vars = [f'{z}_dist_anchor' for z in z_vars]
        avail_dist = [v for v in z_dist_vars if v in udf.columns]
        if avail_dist:
            udf_no_anchor = udf[udf['iso3'] != anchor]
            r = run_panel_gls(udf_no_anchor, 'ca_gdp',
                              avail_dist + controls,
                              f'{union_key}: dist {anchor}')
            if r:
                all_results.append(r)

        # Print demographic profile
        latest = udf[udf['year'] == udf['year'].max()]
        if 'old_dep' in latest.columns and 'youth_dep' in latest.columns:
            print(f"\n  Latest demographics ({int(udf['year'].max())}):")
            print(f"    Old dep: {latest['old_dep'].mean():.1f}% "
                  f"(range {latest['old_dep'].min():.1f}-{latest['old_dep'].max():.1f})")
            print(f"    Youth dep: {latest['youth_dep'].mean():.1f}% "
                  f"(range {latest['youth_dep'].min():.1f}-{latest['youth_dep'].max():.1f})")

    if all_results:
        write_table(all_results, "phase10_within_union_deviations.md",
                    "Within-Union Z Deviations → CA Deviations",
                    note=("*Panel GLS with country and year fixed effects. "
                          "Deviations from union-year cross-sectional means. "
                          "'dist X' = distance from anchor country X.*\n"
                          "*\\*p<0.1, \\*\\*p<0.05, \\*\\*\\*p<0.01*"))

    return all_results


# ── 4. CFA Zone Deep Dive ────────────────────────────────────────────

def cfa_deep_dive(df):
    """Detailed analysis of CFA zone: WAEMU vs CEMAC, projections."""
    print("\n" + "=" * 60)
    print("4. CFA ZONE DEEP DIVE")
    print("=" * 60)

    z_vars = ['Z_1', 'Z_2', 'Z_3']
    controls = ['fiscal_bal_gdp', 'nfa_gdp_lag', 'rgdp_growth']

    cfa_df = get_union_df(df, CFA_JOIN)
    print(f"  CFA sample: {len(cfa_df)} obs, {cfa_df['iso3'].nunique()} countries")

    results = []

    # M1: Full CFA
    r = run_panel_gls(cfa_df, 'ca_gdp', z_vars + controls, 'CFA: Full')
    if r:
        results.append(r)

    # M2: WAEMU only
    waemu_df = get_union_df(df, WAEMU_JOIN)
    r = run_panel_gls(waemu_df, 'ca_gdp', z_vars + controls, 'WAEMU')
    if r:
        results.append(r)

    # M3: CEMAC only
    cemac_df = get_union_df(df, CEMAC_JOIN)
    r = run_panel_gls(cemac_df, 'ca_gdp', z_vars + controls, 'CEMAC')
    if r:
        results.append(r)

    # M4: Age decomposition
    age_vars = ['old_dep', 'youth_dep']
    avail_age = [v for v in age_vars if v in cfa_df.columns]
    if avail_age:
        r = run_panel_gls(cfa_df, 'ca_gdp', avail_age + controls, 'CFA: Age')
        if r:
            results.append(r)

    # M5: WAEMU vs CEMAC interaction
    cfa_df['is_waemu'] = cfa_df['iso3'].isin(WAEMU_ISO3).astype(float)
    cfa_df['Z_1_x_waemu'] = cfa_df['Z_1'] * cfa_df['is_waemu']
    r = run_panel_gls(cfa_df, 'ca_gdp',
                      z_vars + controls + ['is_waemu', 'Z_1_x_waemu'],
                      'CFA: WAEMU int')
    if r:
        results.append(r)

    if results:
        write_table(results, "phase10_cfa_deep_dive.md",
                    "CFA Franc Zone Deep Dive: WAEMU vs CEMAC",
                    note=("*Panel GLS with country and year fixed effects. "
                          "CFA members post-accession. WAEMU: 8 West African "
                          "countries. CEMAC: 6 Central African countries.*\n"
                          "*\\*p<0.1, \\*\\*p<0.05, \\*\\*\\*p<0.01*"))

    # ── CFA Z₁ Projections ──
    print("\n  --- CFA Z₁ Projections ---")
    fp = pd.read_csv(MULTILATERAL_DATA / "full_panel.csv")
    cfa_fp = fp[fp['iso3'].isin(CFA_ISO3)].copy()

    decades = [2000, 2010, 2020, 2030, 2040, 2050, 2060]
    proj_rows = []
    for iso3 in sorted(CFA_ISO3):
        cdata = cfa_fp[cfa_fp['iso3'] == iso3]
        for yr in decades:
            row_data = cdata[cdata['year'] == yr]
            if len(row_data) == 0:
                continue
            proj_rows.append({
                'iso3': iso3,
                'year': yr,
                'Z_1': row_data['Z_1'].values[0],
                'subgroup': 'WAEMU' if iso3 in WAEMU_ISO3 else 'CEMAC',
            })

    proj_df = pd.DataFrame(proj_rows)
    if len(proj_df) > 0:
        # Dispersion by decade
        disp = proj_df.groupby('year')['Z_1'].agg(['std', 'mean', 'count'])
        print("\n  CFA Z₁ Dispersion by Decade:")
        print(disp.to_string(float_format='%.4f'))

        # Deviation from CFA mean
        cfa_mean = proj_df.groupby('year')['Z_1'].mean().rename('Z_1_cfa_mean')
        proj_df = proj_df.merge(cfa_mean, on='year', how='left')
        proj_df['Z_1_dev'] = proj_df['Z_1'] - proj_df['Z_1_cfa_mean']

        # Country × decade pivot
        pivot = proj_df.pivot(index='iso3', columns='year', values='Z_1')
        pivot.columns = [f'Z1_{int(c)}' for c in pivot.columns]
        print("\n  CFA Z₁ Levels by Country:")
        print(pivot.to_string(float_format='%.3f'))

        # Write projection summary
        lines = ["# CFA Franc Zone Z₁ Projections to 2060\n"]
        lines.append("## Cross-Sectional Dispersion\n")
        lines.append("| Decade | Z₁ Mean | Z₁ Std Dev | N Countries |")
        lines.append("|---:|---:|---:|---:|")
        for yr, row in disp.iterrows():
            lines.append(f"| {int(yr)} | {row['mean']:.4f} | "
                         f"{row['std']:.4f} | {int(row['count'])} |")

        lines.append("\n## Country Z₁ Levels\n")
        avail_decades = [d for d in decades if f'Z1_{d}' in pivot.columns]
        lines.append("| Country | Sub | " +
                      " | ".join([str(d) for d in avail_decades]) + " |")
        lines.append("|:---|:---|" +
                      "|".join(["---:" for _ in avail_decades]) + "|")
        for iso3 in sorted(pivot.index):
            sub = 'W' if iso3 in WAEMU_ISO3 else 'C'
            row_str = f"| {iso3} | {sub} |"
            for d in avail_decades:
                col = f'Z1_{d}'
                if col in pivot.columns and pd.notna(pivot.loc[iso3, col]):
                    row_str += f" {pivot.loc[iso3, col]:.3f} |"
                else:
                    row_str += " |"
            lines.append(row_str)

        lines.append("\n*W=WAEMU, C=CEMAC. Z₁ from UN WPP population projections. "
                     "CFA countries are younger (lower Z₁) than EMU — testing "
                     "whether amplification works for youth-driven divergence.*")

        path = OUT_TABLES / "phase10_cfa_projections.md"
        path.write_text('\n'.join(lines))
        print(f"\n  Saved: {path}")

    return results, proj_df if len(proj_df) > 0 else None


# ── 5. Cross-Union Comparison Summary ─────────────────────────────────

def cross_union_summary(df, within_results):
    """Summarize key results across all unions."""
    print("\n" + "=" * 60)
    print("5. CROSS-UNION COMPARISON SUMMARY")
    print("=" * 60)

    z_vars = ['Z_1', 'Z_2', 'Z_3']
    controls = ['fiscal_bal_gdp', 'nfa_gdp_lag', 'rgdp_growth']

    # Compare union members vs non-union floaters
    results = []

    # All union members pooled
    all_union_iso3 = set()
    all_union_join = {}
    for u in UNIONS.values():
        all_union_iso3 |= u['iso3']
        all_union_join.update(u['members'])

    union_df = get_union_df(df, all_union_join)
    r = run_panel_gls(union_df, 'ca_gdp', z_vars + controls,
                      'All unions')
    if r:
        results.append(r)

    # Non-union countries (floaters + non-union peggers)
    non_union = df[~df['iso3'].isin(all_union_iso3)].copy()
    r = run_panel_gls(non_union, 'ca_gdp', z_vars + controls,
                      'Non-union')
    if r:
        results.append(r)

    # Non-union OECD floaters specifically (benchmark)
    oecd_floaters = df[
        (df['iso3'].isin(OECD)) &
        (~df['iso3'].isin(EUROZONE_ISO3)) &
        (df['oecd_floater'] == 1 if 'oecd_floater' in df.columns else True)
    ].copy()
    r = run_panel_gls(oecd_floaters, 'ca_gdp', z_vars + controls,
                      'OECD floaters')
    if r:
        results.append(r)

    if results:
        write_table(results, "phase10_union_vs_nonunion.md",
                    "Union Members vs Non-Union: Demographic CA Effects",
                    note=("*Panel GLS with country and year fixed effects. "
                          "'All unions' pools EMU, CFA, ECCU, CMA members "
                          "post-accession. 'Non-union' excludes all union members.*\n"
                          "*\\*p<0.1, \\*\\*p<0.05, \\*\\*\\*p<0.01*"))

    # Amplification ratio summary
    print("\n  --- Amplification Ratios ---")
    print("  (Union Z₁ coefficient / Non-union Z₁ coefficient)")

    nonunion_z1 = None
    for r in results:
        if r and r['model'] == 'Non-union' and 'Z_1_coef' in r:
            nonunion_z1 = r['Z_1_coef']
            break

    if nonunion_z1 and abs(nonunion_z1) > 0.001:
        for union_key in ['EMU', 'CFA', 'ECCU', 'CMA']:
            if union_key in within_results:
                for wr in within_results[union_key]:
                    if wr and 'Z+ctrl' in wr['model'] and 'Z_1_coef' in wr:
                        ratio = wr['Z_1_coef'] / nonunion_z1
                        sig = stars(wr['Z_1_p'])
                        print(f"  {union_key:6s}: Z₁={wr['Z_1_coef']:8.2f}{sig} "
                              f"/ {nonunion_z1:.2f} = {ratio:.1f}x amplification")
                        break

    # Summary table of Z₁ coefficients across all specifications
    lines = ["# Cross-Union Comparison: Z₁ Coefficient Summary\n"]
    lines.append("| Union | Specification | Z₁ Coef | SE | p-value | N | R² |")
    lines.append("|:---|:---|---:|---:|---:|---:|---:|")

    for union_key in ['EMU', 'CFA', 'ECCU', 'CMA']:
        if union_key not in within_results:
            continue
        for wr in within_results[union_key]:
            if wr and 'Z_1_coef' in wr:
                sig = stars(wr['Z_1_p'])
                lines.append(
                    f"| {union_key} | {wr['model']} | "
                    f"{wr['Z_1_coef']:.2f}{sig} | "
                    f"{wr['Z_1_se']:.2f} | "
                    f"{wr['Z_1_p']:.4f} | "
                    f"{wr['n_obs']} | "
                    f"{wr['r_squared']:.3f} |"
                )

    # Add non-union benchmarks
    for r in results:
        if r and 'Z_1_coef' in r:
            sig = stars(r['Z_1_p'])
            lines.append(
                f"| — | {r['model']} | "
                f"{r['Z_1_coef']:.2f}{sig} | "
                f"{r['Z_1_se']:.2f} | "
                f"{r['Z_1_p']:.4f} | "
                f"{r['n_obs']} | "
                f"{r['r_squared']:.3f} |"
            )

    lines.append("\n*Amplification = union Z₁ coefficient is larger (in absolute "
                 "value) than non-union benchmark, indicating that eliminating "
                 "the exchange rate absorber amplifies demographic imbalances.*")

    path = OUT_TABLES / "phase10_z1_comparison.md"
    path.write_text('\n'.join(lines))
    print(f"\n  Saved: {path}")


# ── 6. Regime Strain for Non-EMU Unions ──────────────────────────────

def non_emu_regime_strain(df, within_results):
    """Project regime strain for CFA, ECCU, CMA using deviation coefficients."""
    print("\n" + "=" * 60)
    print("6. REGIME STRAIN PROJECTIONS FOR NON-EMU UNIONS")
    print("=" * 60)

    fp = pd.read_csv(MULTILATERAL_DATA / "full_panel.csv")
    z_vars = ['Z_1', 'Z_2', 'Z_3']
    decades = [2000, 2010, 2020, 2030, 2040, 2050, 2060]

    for union_key in ['CFA', 'ECCU', 'CMA']:
        union_info = UNIONS[union_key]
        print(f"\n  --- {union_info['label']} ---")

        # Find the dev+ctrl result (or dev-only) for this union
        # from section 3 results
        # We need to re-run deviations to get coefficients
        udf = get_union_df(df, union_info['members'])
        if len(udf) < 30:
            print(f"  Insufficient data ({len(udf)} obs), skipping")
            continue

        u_means = udf.groupby('year')[z_vars + ['ca_gdp']].mean()
        u_means.columns = [f'{c}_union_mean' for c in u_means.columns]
        udf = udf.merge(u_means, on='year', how='left')
        for z in z_vars:
            udf[f'{z}_dev'] = udf[z] - udf[f'{z}_union_mean']
        udf['ca_dev'] = udf['ca_gdp'] - udf['ca_gdp_union_mean']

        z_dev_vars = [f'{z}_dev' for z in z_vars]
        r = run_panel_gls(udf, 'ca_dev', z_dev_vars,
                          f'{union_key}: dev')
        if r is None:
            print(f"  Could not estimate deviation model, skipping strain")
            continue

        z_coefs = {z: r[f'{z}_coef'] for z in z_dev_vars}

        # Get projections
        union_fp = fp[fp['iso3'].isin(union_info['iso3'])].copy()
        proj_rows = []
        for iso3 in sorted(union_info['iso3']):
            cdata = union_fp[union_fp['iso3'] == iso3]
            for yr in decades:
                row_data = cdata[cdata['year'] == yr]
                if len(row_data) == 0:
                    continue
                proj_rows.append({
                    'iso3': iso3, 'year': yr,
                    **{z: row_data[z].values[0] for z in z_vars},
                })

        if not proj_rows:
            print(f"  No projection data for {union_key}")
            continue

        proj_df = pd.DataFrame(proj_rows)

        # Compute union mean and deviations per decade
        for z in z_vars:
            u_decade_mean = proj_df.groupby('year')[z].mean().rename(f'{z}_union_mean')
            proj_df = proj_df.merge(u_decade_mean, on='year', how='left')
            proj_df[f'{z}_dev'] = proj_df[z] - proj_df[f'{z}_union_mean']

        # Compute predicted CA deviation
        strain_rows = []
        for iso3 in sorted(union_info['iso3']):
            for yr in decades:
                cyr = proj_df[(proj_df['iso3'] == iso3) & (proj_df['year'] == yr)]
                if len(cyr) == 0:
                    continue
                predicted_ca_dev = sum(
                    z_coefs[f'{z}_dev'] * cyr[f'{z}_dev'].values[0]
                    for z in z_vars
                )
                strain_rows.append({
                    'iso3': iso3, 'year': yr,
                    'Z_1_dev': cyr['Z_1_dev'].values[0],
                    'predicted_ca_dev': predicted_ca_dev,
                    'strain': abs(predicted_ca_dev),
                })

        strain_df = pd.DataFrame(strain_rows)
        if len(strain_df) == 0:
            continue

        # Display
        target_yr = 2040
        strain_tgt = strain_df[strain_df['year'] == target_yr].sort_values(
            'predicted_ca_dev', ascending=True)
        if len(strain_tgt) > 0:
            print(f"\n  Regime Strain ({target_yr}):")
            for _, row in strain_tgt.iterrows():
                direction = 'surplus' if row['predicted_ca_dev'] > 0 else 'deficit'
                print(f"    {row['iso3']}: CA_dev={row['predicted_ca_dev']:+.2f} "
                      f"({direction})")

        # Write markdown
        pivot = strain_df.pivot(index='iso3', columns='year',
                                values='predicted_ca_dev')
        pivot.columns = [str(int(c)) for c in pivot.columns]

        lines = [f"# {union_info['label']} Regime Strain Index\n"]
        lines.append("Predicted CA/GDP deviation from union mean using Z "
                     "deviation coefficients.\n")
        lines.append("| Country | " + " | ".join(pivot.columns) + " |")
        lines.append("|:---|" + "|".join(["---:" for _ in pivot.columns]) + "|")
        for iso3 in (strain_tgt['iso3'].tolist() if len(strain_tgt) > 0
                     else sorted(pivot.index)):
            if iso3 not in pivot.index:
                continue
            row_str = f"| {iso3} |"
            for col in pivot.columns:
                val = pivot.loc[iso3, col]
                if pd.notna(val):
                    row_str += f" {val:+.2f} |"
                else:
                    row_str += " |"
            lines.append(row_str)

        lines.append(f"\n*Coefficients from within-{union_key} Z deviation "
                     f"regression (N={r['n_obs']}). "
                     f"Z projections from UN WPP.*")

        path = OUT_TABLES / f"phase10_{union_key.lower()}_regime_strain.md"
        path.write_text('\n'.join(lines))
        print(f"  Saved: {path}")


# ── 7. Robustness Checks ─────────────────────────────────────────────

def robustness_checks(df):
    """Comprehensive robustness checks for the monetary union amplification result."""
    print("\n" + "=" * 60)
    print("7. ROBUSTNESS CHECKS")
    print("=" * 60)

    z_vars = ['Z_1', 'Z_2', 'Z_3']
    controls = ['fiscal_bal_gdp', 'nfa_gdp_lag', 'rgdp_growth']

    # ── 7a. CFA excluding oil exporters ──
    print("\n  --- 7a. CFA Excluding Oil Exporters ---")
    # GAB, GNQ, COG, TCD are significant oil/resource exporters
    OIL_EXPORTERS = {'GAB', 'GNQ', 'COG', 'TCD', 'CMR'}
    cfa_no_oil_join = {k: v for k, v in CFA_JOIN.items() if k not in OIL_EXPORTERS}
    cfa_no_oil = get_union_df(df, cfa_no_oil_join)

    oil_results = []

    r = run_panel_gls(cfa_no_oil, 'ca_gdp', z_vars + controls,
                      'CFA excl oil')
    if r:
        oil_results.append(r)

    # CFA oil-only
    cfa_oil_join = {k: v for k, v in CFA_JOIN.items() if k in OIL_EXPORTERS}
    cfa_oil = get_union_df(df, cfa_oil_join)
    r = run_panel_gls(cfa_oil, 'ca_gdp', z_vars + controls,
                      'CFA oil only')
    if r:
        oil_results.append(r)

    # Full CFA for comparison
    cfa_full = get_union_df(df, CFA_JOIN)
    r = run_panel_gls(cfa_full, 'ca_gdp', z_vars + controls, 'CFA full')
    if r:
        oil_results.append(r)

    if oil_results:
        write_table(oil_results, "phase10_robustness_cfa_oil.md",
                    "Robustness: CFA Excluding Oil Exporters",
                    note=("*Panel GLS with country and year fixed effects. "
                          "Oil exporters: GAB, GNQ, COG, TCD, CMR.*\n"
                          "*\\*p<0.1, \\*\\*p<0.05, \\*\\*\\*p<0.01*"))

    # ── 7b. Time period stability ──
    print("\n  --- 7b. Time Period Stability ---")
    period_results = []
    periods = {
        'Pre-2000': (1960, 1999),
        '2000-2012': (2000, 2012),
        'Post-2012': (2013, 2024),
    }

    for union_key in ['EMU', 'CFA']:
        union_info = UNIONS[union_key]
        udf = get_union_df(df, union_info['members'])
        for period_name, (yr_start, yr_end) in periods.items():
            sub = udf[(udf['year'] >= yr_start) & (udf['year'] <= yr_end)]
            r = run_panel_gls(sub, 'ca_gdp', z_vars + controls,
                              f'{union_key}: {period_name}')
            if r:
                period_results.append(r)

    if period_results:
        write_table(period_results, "phase10_robustness_periods.md",
                    "Robustness: Time Period Stability",
                    note=("*Panel GLS with country and year fixed effects. "
                          "Union members post-accession, split by time period.*\n"
                          "*\\*p<0.1, \\*\\*p<0.05, \\*\\*\\*p<0.01*"))

    # ── 7c. Controlling for trade openness and KAOPEN ──
    print("\n  --- 7c. Additional Controls ---")
    addl_results = []
    extra_controls_sets = {
        'Baseline': controls,
        '+ trade_openness': controls + ['trade_openness'],
        '+ kaopen': controls + ['kaopen'],
        '+ trade + kaopen': controls + ['trade_openness', 'kaopen'],
    }

    for union_key in ['EMU', 'CFA']:
        union_info = UNIONS[union_key]
        udf = get_union_df(df, union_info['members'])
        for ctrl_name, ctrl_vars in extra_controls_sets.items():
            avail = [v for v in ctrl_vars if v in udf.columns]
            r = run_panel_gls(udf, 'ca_gdp', z_vars + avail,
                              f'{union_key}: {ctrl_name}')
            if r:
                addl_results.append(r)

    if addl_results:
        write_table(addl_results, "phase10_robustness_controls.md",
                    "Robustness: Additional Controls",
                    note=("*Panel GLS with country and year fixed effects. "
                          "Tests whether amplification survives inclusion of "
                          "trade openness and capital account openness (KAOPEN).*\n"
                          "*\\*p<0.1, \\*\\*p<0.05, \\*\\*\\*p<0.01*"))

    # ── 7d. Excluding outlier countries ──
    print("\n  --- 7d. Excluding Outlier Countries ---")
    outlier_results = []

    # EMU: exclude GRC, IRL (known CA outliers)
    emu_join_excl = {k: v for k, v in EUROZONE_JOIN.items()
                     if k not in {'GRC', 'IRL'}}
    emu_excl = get_union_df(df, emu_join_excl)
    r = run_panel_gls(emu_excl, 'ca_gdp', z_vars + controls,
                      'EMU excl GRC,IRL')
    if r:
        outlier_results.append(r)

    # EMU: exclude small countries (< 5M pop)
    SMALL_EMU = {'CYP', 'EST', 'LVA', 'LTU', 'LUX', 'MLT', 'SVN', 'SVK'}
    emu_large_join = {k: v for k, v in EUROZONE_JOIN.items()
                      if k not in SMALL_EMU}
    emu_large = get_union_df(df, emu_large_join)
    r = run_panel_gls(emu_large, 'ca_gdp', z_vars + controls,
                      'EMU large only')
    if r:
        outlier_results.append(r)

    # CFA: exclude GNQ (tiny population, massive oil GDP)
    cfa_no_gnq_join = {k: v for k, v in CFA_JOIN.items() if k != 'GNQ'}
    cfa_no_gnq = get_union_df(df, cfa_no_gnq_join)
    r = run_panel_gls(cfa_no_gnq, 'ca_gdp', z_vars + controls,
                      'CFA excl GNQ')
    if r:
        outlier_results.append(r)

    # CFA: WAEMU only (non-oil, homogeneous)
    waemu_df = get_union_df(df, WAEMU_JOIN)
    r = run_panel_gls(waemu_df, 'ca_gdp', z_vars + controls, 'WAEMU only')
    if r:
        outlier_results.append(r)

    if outlier_results:
        write_table(outlier_results, "phase10_robustness_outliers.md",
                    "Robustness: Excluding Outlier Countries",
                    note=("*Panel GLS with country and year fixed effects. "
                          "Tests sensitivity to exclusion of known outliers.*\n"
                          "*\\*p<0.1, \\*\\*p<0.05, \\*\\*\\*p<0.01*"))

    # ── 7e. Non-union peggers as placebo ──
    print("\n  --- 7e. Placebo: Non-Union Peggers ---")
    placebo_results = []

    # Countries that peg but are NOT in a monetary union
    all_union_iso3 = set()
    for u in UNIONS.values():
        all_union_iso3 |= u['iso3']

    # Non-union peggers (regime_3cat == 1 or is_peg == 1)
    if 'is_peg' in df.columns:
        peg_col = 'is_peg'
    elif 'regime_3cat' in df.columns:
        df['is_peg_temp'] = (df['regime_3cat'] == 1).astype(float)
        peg_col = 'is_peg_temp'
    else:
        peg_col = None

    if peg_col:
        non_union_peggers = df[
            (~df['iso3'].isin(all_union_iso3)) &
            (df[peg_col] == 1)
        ].copy()
        r = run_panel_gls(non_union_peggers, 'ca_gdp', z_vars + controls,
                          'Non-union peggers')
        if r:
            placebo_results.append(r)

        # Non-union floaters
        non_union_floaters = df[
            (~df['iso3'].isin(all_union_iso3)) &
            (df[peg_col] == 0)
        ].copy()
        r = run_panel_gls(non_union_floaters, 'ca_gdp', z_vars + controls,
                          'Non-union floaters')
        if r:
            placebo_results.append(r)

    # Full union sample for comparison
    all_union_join = {}
    for u in UNIONS.values():
        all_union_join.update(u['members'])
    union_df = get_union_df(df, all_union_join)
    r = run_panel_gls(union_df, 'ca_gdp', z_vars + controls,
                      'All union members')
    if r:
        placebo_results.append(r)

    if placebo_results:
        write_table(placebo_results, "phase10_robustness_placebo.md",
                    "Placebo: Non-Union Peggers vs Union Members vs Floaters",
                    note=("*Panel GLS with country and year fixed effects. "
                          "Non-union peggers share exchange rate fixity but NOT "
                          "monetary union constraints. If amplification is about "
                          "removing the exchange rate absorber within a union "
                          "(not just pegging), union members should show larger "
                          "demographic effects than non-union peggers.*\n"
                          "*\\*p<0.1, \\*\\*p<0.05, \\*\\*\\*p<0.01*"))

    # ── 7f. Z₁ × KAOPEN within unions ──
    print("\n  --- 7f. Z₁ × KAOPEN Within Unions ---")
    kaopen_results = []

    for union_key in ['EMU', 'CFA', 'CMA']:
        union_info = UNIONS[union_key]
        udf = get_union_df(df, union_info['members'])
        if 'kaopen' not in udf.columns:
            continue
        udf['Z_1_x_kaopen'] = udf['Z_1'] * udf['kaopen']
        r = run_panel_gls(udf, 'ca_gdp',
                          z_vars + controls + ['kaopen', 'Z_1_x_kaopen'],
                          f'{union_key}: Z₁×KA')
        if r:
            kaopen_results.append(r)

    if kaopen_results:
        write_table(kaopen_results, "phase10_robustness_kaopen.md",
                    "Robustness: Z₁ × KAOPEN Within Unions",
                    note=("*Panel GLS with country and year fixed effects. "
                          "Tests whether capital account openness modifies the "
                          "demographic effect within monetary unions.*\n"
                          "*\\*p<0.1, \\*\\*p<0.05, \\*\\*\\*p<0.01*"))

    # ── 7g. Age decomposition across unions ──
    print("\n  --- 7g. Age Decomposition Comparison ---")
    age_results = []
    age_vars = ['old_dep', 'youth_dep']

    for union_key in ['EMU', 'CFA', 'ECCU', 'CMA']:
        union_info = UNIONS[union_key]
        udf = get_union_df(df, union_info['members'])
        avail = [v for v in age_vars if v in udf.columns]
        if avail:
            r = run_panel_gls(udf, 'ca_gdp', avail + controls,
                              f'{union_key}: Age')
            if r:
                age_results.append(r)

    # Non-union benchmark
    all_union_iso3 = set()
    for u in UNIONS.values():
        all_union_iso3 |= u['iso3']
    non_union = df[~df['iso3'].isin(all_union_iso3)]
    avail = [v for v in age_vars if v in non_union.columns]
    if avail:
        r = run_panel_gls(non_union, 'ca_gdp', avail + controls,
                          'Non-union')
        if r:
            age_results.append(r)

    if age_results:
        write_table(age_results, "phase10_robustness_age_decomp.md",
                    "Robustness: Age Decomposition Across Unions",
                    note=("*Panel GLS with country and year fixed effects. "
                          "Uses old-age and youth dependency ratios instead of "
                          "Z demographic composites. Tests whether amplification "
                          "is driven by old-age (EMU) vs youth (CFA) demographics.*\n"
                          "*\\*p<0.1, \\*\\*p<0.05, \\*\\*\\*p<0.01*"))

    # ── 7h. Winsorized CA ──
    print("\n  --- 7h. Winsorized CA ---")
    wins_results = []

    for union_key in ['EMU', 'CFA']:
        union_info = UNIONS[union_key]
        udf = get_union_df(df, union_info['members']).copy()
        # Winsorize CA at p1/p99
        p1, p99 = udf['ca_gdp'].quantile([0.01, 0.99])
        udf['ca_gdp_wins'] = udf['ca_gdp'].clip(p1, p99)

        r = run_panel_gls(udf, 'ca_gdp_wins', z_vars + controls,
                          f'{union_key}: CA wins')
        if r:
            wins_results.append(r)

    if wins_results:
        write_table(wins_results, "phase10_robustness_winsorized.md",
                    "Robustness: Winsorized CA/GDP (p1/p99)",
                    note=("*Panel GLS with country and year fixed effects. "
                          "CA/GDP winsorized at 1st/99th percentile within each "
                          "union to reduce outlier influence.*\n"
                          "*\\*p<0.1, \\*\\*p<0.05, \\*\\*\\*p<0.01*"))

    # ── 7i. Commodity price controls for CFA ──
    print("\n  --- 7i. GDP Per Capita Control (CFA) ---")
    gdp_results = []
    cfa_full = get_union_df(df, CFA_JOIN)

    # GDP per capita as control (proxy for resource-driven income)
    if 'gdp_pc_ppp' in cfa_full.columns:
        cfa_full['log_gdp_pc'] = np.log(cfa_full['gdp_pc_ppp'].clip(lower=100))
        r = run_panel_gls(cfa_full, 'ca_gdp',
                          z_vars + controls + ['log_gdp_pc'],
                          'CFA + log GDP/cap')
        if r:
            gdp_results.append(r)

    # CFA with trade openness control
    if 'trade_openness' in cfa_full.columns:
        r = run_panel_gls(cfa_full, 'ca_gdp',
                          z_vars + controls + ['trade_openness'],
                          'CFA + trade open')
        if r:
            gdp_results.append(r)

    # Both
    extra = []
    if 'gdp_pc_ppp' in cfa_full.columns:
        extra.append('log_gdp_pc')
    if 'trade_openness' in cfa_full.columns:
        extra.append('trade_openness')
    if extra:
        r = run_panel_gls(cfa_full, 'ca_gdp', z_vars + controls + extra,
                          'CFA + all extras')
        if r:
            gdp_results.append(r)

    if gdp_results:
        write_table(gdp_results, "phase10_robustness_cfa_controls.md",
                    "Robustness: CFA with GDP/Capita and Trade Openness Controls",
                    note=("*Panel GLS with country and year fixed effects. "
                          "Tests whether CFA amplification survives controlling "
                          "for income level (resource rents proxy) and trade "
                          "openness.*\n"
                          "*\\*p<0.1, \\*\\*p<0.05, \\*\\*\\*p<0.01*"))

    # ── Summary ──
    print("\n  --- Robustness Summary ---")
    print("  Tables saved:")
    for f in sorted(OUT_TABLES.glob("phase10_robustness_*.md")):
        print(f"    {f.name}")


# ── Main ──────────────────────────────────────────────────────────────

def main():
    print("=" * 70)
    print("PHASE 10: MONETARY UNION COMPARISON — BEYOND THE EUROZONE")
    print("=" * 70)

    df = pd.read_csv(DATA / "trilemma_panel.csv")
    print(f"Panel: {len(df)} obs, {df['iso3'].nunique()} countries")
    print(f"Year range: {df['year'].min()}-{df['year'].max()}")

    # Print union coverage
    print("\n  Union coverage in panel:")
    for union_key, union_info in UNIONS.items():
        udf = get_union_df(df, union_info['members'])
        joint = udf.dropna(subset=['ca_gdp', 'Z_1'])
        print(f"    {union_key:6s} ({union_info['label']:30s}): "
              f"{joint['iso3'].nunique():2d} countries, {len(joint):5d} obs")

    # 1. Within-union CA regressions
    within_results = within_union_ca(df)

    # 2. Pooled monetary union analysis
    pdf, pooled_results = pooled_union_analysis(df)

    # 3. Within-union Z deviations
    dev_results = within_union_deviations(df)

    # 4. CFA zone deep dive
    cfa_results, cfa_proj = cfa_deep_dive(df)

    # 5. Cross-union comparison summary
    cross_union_summary(df, within_results)

    # 6. Regime strain for non-EMU unions
    non_emu_regime_strain(df, within_results)

    # 7. Robustness checks
    robustness_checks(df)

    print("\n" + "=" * 70)
    print("PHASE 10 COMPLETE")
    print("=" * 70)


if __name__ == '__main__':
    main()
